###########################
### PANDEMICS AND POLITICAL DEVELOPMENT
### THE ELECTORAL LEGACY OF THE BLACK DEATH IN GERMANY
### WORLD POLITICS, VOL. 73, NO. 3 (JULY 2021)
### REPLICATION FILE OF THE EMPIRICAL ANALYSIS
### DANIEL W. GINGERICH & JAN P. VOGLER
###########################

###########################
### REPLICATION FILES PART 4
### EARLY NINETEENTH-CENTURY PRUSSIA
###########################

###########################
### INSTALL REQUIRED PACKAGES
###########################

install.packages("geosphere")
install.packages("AER")



###########################
### DATA MANAGEMENT
###########################

###########################
### DATA ON THE NUMBER OF FARMS/ESTATES OF DIFFERENT SIZES
### BY BECKER ET AL. (2014)
### IFO PRUSSIAN ECONOMIC HISTORY DATABASE (IPEHD)
###########################

data_outcomes_01=read.csv("ipehd_1816_agri_land.csv")

# AGGREGATE DIFFERENT MEASURES (SUM) BY KREISKEY1800 TO OBTAIN VALUES FOR KK1800 UNITS

data_outcomes_01_agg_ov300=aggregate(formula = land1816_far_ext_ov300 ~ kreiskey1800, data = data_outcomes_01, FUN = sum) # FARMS WITH MORE THAN 300 PRUSSIAN MORGEN
data_outcomes_01_agg_15to30=aggregate(formula = land1816_far_ext_15to30 ~ kreiskey1800, data = data_outcomes_01, FUN = sum) # FARMS WITH 15 TO 300 PRUSSIAN MORGEN
data_outcomes_01_agg_un15=aggregate(formula = land1816_far_ext_un15 ~ kreiskey1800, data = data_outcomes_01, FUN = sum) # FARMS WITH LESS THAN 15 PRUSSIAN MORGEN



###########################
### DATA ON POPULATION DEMOGRAPHICS IN 1816
### BY BECKER ET AL. (2014)
### IFO PRUSSIAN ECONOMIC HISTORY DATABASE (IPEHD)
###########################

data_outcomes_02=read.csv("ipehd_1816_pop_demo.csv")

data_outcomes_02=data_outcomes_02[,c("kreiskey1800","county","rb","id1816","pop1816_tot")] # NEED TOTAL POPULATION

# AGGREGATE TOTAL POPULATION (SUM) BY KREISKEY1800 TO OBTAIN VALUES FOR KK1800 UNITS

data_outcomes_02_agg_totpop=aggregate(formula = pop1816_tot ~ kreiskey1800, data = data_outcomes_02, FUN = sum)



###########################
### DATA ON THE NUMBER OF AGRICULTURAL SERVANTS IN 1819
### BY BECKER ET AL. (2014)
### IFO PRUSSIAN ECONOMIC HISTORY DATABASE (IPEHD)
###########################

data_outcomes_03=read.csv("ipehd_1819_occ_agri.csv")

data_outcomes_03=data_outcomes_03[,c("kreiskey1800","county","rb","id1819","occ1819_servant_m_farm","occ1819_servant_f_farm")]

# AGGREGATE DIFFERENT MEASURES (SUM) BY KREISKEY1800 TO OBTAIN VALUES FOR KK1800 UNITS

data_outcomes_03_agg_mservant=aggregate(formula = occ1819_servant_m_farm ~ kreiskey1800, data = data_outcomes_03, FUN = sum)
data_outcomes_03_agg_fservant=aggregate(formula = occ1819_servant_f_farm ~ kreiskey1800, data = data_outcomes_03, FUN = sum)



###########################
### DATA ON CONTROL VARIABLES
###########################

data_control_var=read.csv("data_pru_control_var.csv")



###########################
### DATA ON THE GEOGRAPHIC LOCATION OF PRUSSIAN COUNTIES
###########################

data_county_geoloc=read.csv("data_pru_location.csv")



###########################
### DATA ON MERGING PROPERTIES
### CONTAINING INFORMATION ON HOW TO ASSIGN KK1871 UNITS TO KK1800 UNITS
###########################

data_merging_info=read.csv("ipehd_merge_county.csv")

length(unique(data_merging_info$kreiskey1800)) # 281 UNIQUE KREISKEY1800 VALUES
summary(data_merging_info$kreiskey1800) # RANGE: 1-281

length(unique(data_merging_info$kreiskey1871)) # 453 UNIQUE KREISKEY1871 VALUES (1 NA)
summary(data_merging_info$kreiskey1871) # RANGE 1-453



###########################
### CREATE 1871 MERGING FILE
### CONTAINING 1 OBSERVATION PER KREISKEY1871 UNIT
### INFORMATION ON ASSOCIATED KREISKEY1800 UNIT
###########################

data_merging_1871=data_merging_info[,c("kreiskey1800","kreiskey1871","county1871")]

data_merging_1871=data_merging_1871[-which(is.na(data_merging_1871$kreiskey1871)),]

data_merging_1871=data_merging_1871[-which(is.na(data_merging_1871$kreiskey1800)),]

data_merging_1871=unique(data_merging_1871)

data_merging_1871=data_merging_1871[-which(data_merging_1871$kreiskey1871==224)[2],] # KREISKEY1871 IS ASSIGNED TO TWO DIFFERENT KREISKEY1800 UNITS; IT IS NECESSARY LIMIT IT TO ONE



###########################
### MERGE WITH LOCATION DATA
###########################

data_location_1871=merge(data_merging_1871,data_county_geoloc, by=c("kreiskey1871","kreiskey1800","county1871"))

# IDENTIFY AN APPROXIMATE LATITUDE AND LONGITUDE VALUE FOR EACH AGGREGATE KREISKEY1800 UNIT
# NOTE: AGGREGATION BY THE MEAN ALLOWS US TO IDENTIFY ONLY AN APPROXIMATE (AND IMPERFECT) VALUE FOR BOTH LATITUDE AND LONGITUDE
# CONSIDERING THAT THE KREISKEY1800 UNITS ARE BASED ON UNITS IN VERY CLOSE GEOGRAPHIC PROXIMITY, THIS IMPERFECT APPROXIMATION WILL BE SUFFICIENT FOR OUR PURPOSES

data_location_agg_lat=aggregate(formula = Lat ~ kreiskey1800, data = data_location_1871, FUN = mean)
data_location_agg_long=aggregate(formula = Long ~ kreiskey1800, data = data_location_1871, FUN = mean)



###########################
### DATA ON PROXIMITY TO GEOGRAPHIC FEATURES
###########################

data_geography=read.csv("data_pru_geography.csv")

data_geography=data_geography[,c("kreiskey1871","dist_ocean_km","dist_majport_km","dist_river_km","dist_tradecity_km")]



###########################
### MERGE WITH GEOGRAPHIC DATA
###########################

data_geography_1871=merge(data_merging_1871,data_geography, by=c("kreiskey1871"))

# IDENTIFY AVERAGE DISTANCES TO GEOGRAPHIC FEATURES FOR EACH KREISKEY1800 UNIT

data_geography_agg_ocean=aggregate(formula = dist_ocean_km ~ kreiskey1800, data = data_geography_1871, FUN = mean)
data_geography_agg_port=aggregate(formula = dist_majport_km ~ kreiskey1800, data = data_geography_1871, FUN = mean)
data_geography_agg_river=aggregate(formula = dist_river_km ~ kreiskey1800, data = data_geography_1871, FUN = mean)
data_geography_agg_trade=aggregate(formula = dist_tradecity_km ~ kreiskey1800, data = data_geography_1871, FUN = mean)



###########################
### CREATE FINAL DATASET FOR THE EMPIRICAL ANALYSIS
###########################

kreiskey1800=1:281

data_prussia_fin=as.data.frame(x=kreiskey1800,)

data_prussia_fin=merge(data_prussia_fin, data_outcomes_01_agg_ov300, by="kreiskey1800", all.x=T)
data_prussia_fin=merge(data_prussia_fin, data_outcomes_01_agg_15to30, by="kreiskey1800", all.x=T)
data_prussia_fin=merge(data_prussia_fin, data_outcomes_01_agg_un15, by="kreiskey1800", all.x=T)
data_prussia_fin=merge(data_prussia_fin, data_outcomes_02_agg_totpop, by="kreiskey1800", all.x=T)
data_prussia_fin=merge(data_prussia_fin, data_outcomes_03_agg_fservant, by="kreiskey1800", all.x=T)
data_prussia_fin=merge(data_prussia_fin, data_outcomes_03_agg_mservant, by="kreiskey1800", all.x=T)
data_prussia_fin=merge(data_prussia_fin, data_location_agg_lat, by="kreiskey1800", all.x=T)
data_prussia_fin=merge(data_prussia_fin, data_location_agg_long, by="kreiskey1800", all.x=T)
data_prussia_fin=merge(data_prussia_fin, data_geography_agg_ocean, by="kreiskey1800", all.x=T)
data_prussia_fin=merge(data_prussia_fin, data_geography_agg_river, by="kreiskey1800", all.x=T)
data_prussia_fin=merge(data_prussia_fin, data_geography_agg_port, by="kreiskey1800", all.x=T)
data_prussia_fin=merge(data_prussia_fin, data_geography_agg_trade, by="kreiskey1800", all.x=T)



###########################
### CREATE FINAL OUTCOME VARIABLES
###########################

data_prussia_fin$prop_servants_farm=(data_prussia_fin$occ1819_servant_m_farm+data_prussia_fin$occ1819_servant_f_farm)/data_prussia_fin$pop1816_tot

data_prussia_fin$prop_large_farms=data_prussia_fin$land1816_far_ext_ov300/(data_prussia_fin$land1816_far_ext_ov300+data_prussia_fin$land1816_far_ext_15to30+data_prussia_fin$land1816_far_ext_un15)



###########################
### MERGE WITH CONTROL VARIABLES
###########################

data_prussia_fin=merge(data_prussia_fin,data_control_var, by="kreiskey1800")



###########################
### DATA ON BLACK DEATH OUTBREAKS
### BY JEDWAB ET AL. (2019)
###########################

data_bd_outbreaks=read.csv("data_pru_bd_outbreaks.csv", stringsAsFactors=F)



###########################
### CALCULATING BDEI SCORES
###########################

###########################
### COMPUTE THE GLOBAL MAXIMUM DISTANCE
###########################

library(geosphere)
max_store=rep(NA,length(data_prussia_fin$kreiskey1800))
for (j in 1:length(data_prussia_fin$kreiskey1800)){
  store=rep(NA,length(data_bd_outbreaks$int))
  for (i in 1:length(data_bd_outbreaks$int)){
    store[i]=distm(c(data_prussia_fin$Long[j], data_prussia_fin$Lat[j]), c(data_bd_outbreaks$longitude[i], data_bd_outbreaks$latitude[i]), fun = distHaversine)
    max_store[j]=max(store)
  }
}
max_dist=max(max_store)



###########################
### COMPUTE BLACK DEATH EXPOSURE INTENSITY SCORE (V1-5)
###########################

for (j in 1:length(data_prussia_fin$kreiskey1800)){
  store=rep(NA,length(data_bd_outbreaks$int))
  for (i in 1:length(data_bd_outbreaks$int)){
    store[i]=distm(c(data_prussia_fin$Long[j], data_prussia_fin$Lat[j]), c(data_bd_outbreaks$longitude[i], data_bd_outbreaks$latitude[i]), fun = distHaversine)
    data_prussia_fin$BDEI_score[j]=sum(data_bd_outbreaks$int*(1-store/max_dist)^3)
    data_prussia_fin$BDEI_score2[j]=sum(data_bd_outbreaks$int*(1-store/max_dist)^6)
    data_prussia_fin$BDEI_score3[j]=sum(data_bd_outbreaks$int*(1-store/max_dist)^9)
    data_prussia_fin$BDEI_score4[j]=sum(data_bd_outbreaks$int*(1-store/max_dist)^12)
    data_prussia_fin$BDEI_score5[j]=sum(data_bd_outbreaks$int*(1-store/max_dist)^15)
  }
}



###########################
### STANDARDIZE BLACK DEATH EXPOSURE INTENSITY SCORE (V1-5)
###########################

data_prussia_fin$BDEI_score=scale(data_prussia_fin$BDEI_score)[,1]
data_prussia_fin$BDEI_score2=scale(data_prussia_fin$BDEI_score2)[,1]
data_prussia_fin$BDEI_score3=scale(data_prussia_fin$BDEI_score3)[,1]
data_prussia_fin$BDEI_score4=scale(data_prussia_fin$BDEI_score4)[,1]
data_prussia_fin$BDEI_score5=scale(data_prussia_fin$BDEI_score5)[,1]



###########################
### MAIN EMPIRICAL ANALYSIS IN THE ARTICLE
###########################

###########################
### PROPORTION OF LARGE ESTATES (1816)
### TABLE 8 IN THE ARTICLE
###########################

prop_large_farms_lm1=lm(prop_large_farms ~ BDEI_score, data=data_prussia_fin)
summary(prop_large_farms_lm1)

prop_large_farms_lm2=lm(prop_large_farms ~ BDEI_score2, data=data_prussia_fin)
summary(prop_large_farms_lm2)

prop_large_farms_lm3=lm(prop_large_farms ~ BDEI_score3, data=data_prussia_fin)
summary(prop_large_farms_lm3)

prop_large_farms_lm4=lm(prop_large_farms ~ BDEI_score4, data=data_prussia_fin)
summary(prop_large_farms_lm4)

prop_large_farms_lm5=lm(prop_large_farms ~ BDEI_score5, data=data_prussia_fin)
summary(prop_large_farms_lm5)

prop_large_farms_ctrl_lm1=lm(prop_large_farms ~ BDEI_score + urban_density_norm + dist_majport_km + dist_tradecity_km + dist_ocean_km + dist_river_km + elevation, data=data_prussia_fin)
summary(prop_large_farms_ctrl_lm1)

prop_large_farms_ctrl_lm2=lm(prop_large_farms ~ BDEI_score2 + urban_density_norm + dist_majport_km + dist_tradecity_km + dist_ocean_km + dist_river_km + elevation, data=data_prussia_fin)
summary(prop_large_farms_ctrl_lm2)

prop_large_farms_ctrl_lm3=lm(prop_large_farms ~ BDEI_score3 + urban_density_norm + dist_majport_km + dist_tradecity_km + dist_ocean_km + dist_river_km + elevation, data=data_prussia_fin)
summary(prop_large_farms_ctrl_lm3)

prop_large_farms_ctrl_lm4=lm(prop_large_farms ~ BDEI_score4 + urban_density_norm + dist_majport_km + dist_tradecity_km + dist_ocean_km + dist_river_km + elevation, data=data_prussia_fin)
summary(prop_large_farms_ctrl_lm4)

prop_large_farms_ctrl_lm5=lm(prop_large_farms ~ BDEI_score5 + urban_density_norm + dist_majport_km + dist_tradecity_km + dist_ocean_km + dist_river_km + elevation, data=data_prussia_fin)
summary(prop_large_farms_ctrl_lm5)



###########################
### PROPORTION OF AGRICULTURAL SERVANTS (OF TOTAL POPULATION) (1816/1819)
### TABLE 9 IN THE ARTICLE
###########################

servants_lm1=lm(prop_servants_farm ~ BDEI_score, data=data_prussia_fin)
summary(servants_lm1)

servants_lm2=lm(prop_servants_farm ~ BDEI_score2, data=data_prussia_fin)
summary(servants_lm2)

servants_lm3=lm(prop_servants_farm ~ BDEI_score3, data=data_prussia_fin)
summary(servants_lm3)

servants_lm4=lm(prop_servants_farm ~ BDEI_score4, data=data_prussia_fin)
summary(servants_lm4)

servants_lm5=lm(prop_servants_farm ~ BDEI_score5, data=data_prussia_fin)
summary(servants_lm5)

servants_ctrl_lm1=lm(prop_servants_farm ~ BDEI_score + urban_density_norm + dist_majport_km + dist_tradecity_km + dist_ocean_km + dist_river_km + elevation, data=data_prussia_fin)
summary(servants_ctrl_lm1)

servants_ctrl_lm2=lm(prop_servants_farm ~ BDEI_score2 + urban_density_norm + dist_majport_km + dist_tradecity_km + dist_ocean_km + dist_river_km + elevation, data=data_prussia_fin)
summary(servants_ctrl_lm2)

servants_ctrl_lm3=lm(prop_servants_farm ~ BDEI_score3 + urban_density_norm + dist_majport_km + dist_tradecity_km + dist_ocean_km + dist_river_km + elevation, data=data_prussia_fin)
summary(servants_ctrl_lm3)

servants_ctrl_lm4=lm(prop_servants_farm ~ BDEI_score4 + urban_density_norm + dist_majport_km + dist_tradecity_km + dist_ocean_km + dist_river_km + elevation, data=data_prussia_fin)
summary(servants_ctrl_lm4)

servants_ctrl_lm5=lm(prop_servants_farm ~ BDEI_score5 + urban_density_norm + dist_majport_km + dist_tradecity_km + dist_ocean_km + dist_river_km + elevation, data=data_prussia_fin)
summary(servants_ctrl_lm5)



###########################
### ADDITIONAL EMPIRICAL ANALYSIS IN THE SUPPLEMENTARY MATERIAL
###########################

###########################
### DESCRIPTIVE SUMMARY STATISTICS
### TABLE A29 IN THE SUPPLEMENTARY MATERIAL
###########################

summary(data_prussia_fin)



###########################
### TOBIT MODELS AS AN ALTERNATIVE SPECIFICATION
###########################

###########################
### PROPORTION OF LARGE ESTATES (1816)
### TABLE A30 IN THE SUPPLEMENTARY MATERIAL
###########################

library(AER)

prop_large_farms_tobit1=tobit(prop_large_farms ~ BDEI_score, left=0, right=1, data=data_prussia_fin)
summary(prop_large_farms_tobit1)

prop_large_farms_tobit2=tobit(prop_large_farms ~ BDEI_score2, left=0, right=1, data=data_prussia_fin)
summary(prop_large_farms_tobit2)

prop_large_farms_tobit3=tobit(prop_large_farms ~ BDEI_score3, left=0, right=1, data=data_prussia_fin)
summary(prop_large_farms_tobit3)

prop_large_farms_tobit4=tobit(prop_large_farms ~ BDEI_score4, left=0, right=1, data=data_prussia_fin)
summary(prop_large_farms_tobit4)

prop_large_farms_tobit5=tobit(prop_large_farms ~ BDEI_score5, left=0, right=1, data=data_prussia_fin)
summary(prop_large_farms_tobit5)

prop_large_farms_ctrl_tobit1=tobit(prop_large_farms ~ BDEI_score + urban_density_norm + dist_majport_km + dist_tradecity_km + dist_ocean_km + dist_river_km + elevation, left=0, right=1, data=data_prussia_fin)
summary(prop_large_farms_ctrl_tobit1)

prop_large_farms_ctrl_tobit2=tobit(prop_large_farms ~ BDEI_score2 + urban_density_norm + dist_majport_km + dist_tradecity_km + dist_ocean_km + dist_river_km + elevation, left=0, right=1, data=data_prussia_fin)
summary(prop_large_farms_ctrl_tobit2)

prop_large_farms_ctrl_tobit3=tobit(prop_large_farms ~ BDEI_score3 + urban_density_norm + dist_majport_km + dist_tradecity_km + dist_ocean_km + dist_river_km + elevation, left=0, right=1, data=data_prussia_fin)
summary(prop_large_farms_ctrl_tobit3)

prop_large_farms_ctrl_tobit4=tobit(prop_large_farms ~ BDEI_score4 + urban_density_norm + dist_majport_km + dist_tradecity_km + dist_ocean_km + dist_river_km + elevation, left=0, right=1, data=data_prussia_fin)
summary(prop_large_farms_ctrl_tobit4)

prop_large_farms_ctrl_tobit5=tobit(prop_large_farms ~ BDEI_score5 + urban_density_norm + dist_majport_km + dist_tradecity_km + dist_ocean_km + dist_river_km + elevation, left=0, right=1, data=data_prussia_fin)
summary(prop_large_farms_ctrl_tobit5)



###########################
### PROPORTION OF AGRICULTURAL SERVANTS (OF TOTAL POPULATION) (1816/1819)
### TABLE A31 IN THE SUPPLEMENTARY MATERIAL
###########################

servants_tobit1=tobit(prop_servants_farm ~ BDEI_score, left=0, right=1, data=data_prussia_fin)
summary(servants_tobit1)

servants_tobit2=tobit(prop_servants_farm ~ BDEI_score2, left=0, right=1, data=data_prussia_fin)
summary(servants_tobit2)

servants_tobit3=tobit(prop_servants_farm ~ BDEI_score3, left=0, right=1, data=data_prussia_fin)
summary(servants_tobit3)

servants_tobit4=tobit(prop_servants_farm ~ BDEI_score4, left=0, right=1, data=data_prussia_fin)
summary(servants_tobit4)

servants_tobit5=tobit(prop_servants_farm ~ BDEI_score5, left=0, right=1, data=data_prussia_fin)
summary(servants_tobit5)

servants_ctrl_tobit1=tobit(prop_servants_farm ~ BDEI_score + urban_density_norm + dist_majport_km + dist_tradecity_km + dist_ocean_km + dist_river_km + elevation, left=0, right=1, data=data_prussia_fin)
summary(servants_ctrl_tobit1)

servants_ctrl_tobit2=tobit(prop_servants_farm ~ BDEI_score2 + urban_density_norm + dist_majport_km + dist_tradecity_km + dist_ocean_km + dist_river_km + elevation, left=0, right=1, data=data_prussia_fin)
summary(servants_ctrl_tobit2)

servants_ctrl_tobit3=tobit(prop_servants_farm ~ BDEI_score3 + urban_density_norm + dist_majport_km + dist_tradecity_km + dist_ocean_km + dist_river_km + elevation, left=0, right=1, data=data_prussia_fin)
summary(servants_ctrl_tobit3)

servants_ctrl_tobit4=tobit(prop_servants_farm ~ BDEI_score4 + urban_density_norm + dist_majport_km + dist_tradecity_km + dist_ocean_km + dist_river_km + elevation, left=0, right=1, data=data_prussia_fin)
summary(servants_ctrl_tobit4)

servants_ctrl_tobit5=tobit(prop_servants_farm ~ BDEI_score5 + urban_density_norm + dist_majport_km + dist_tradecity_km + dist_ocean_km + dist_river_km + elevation, left=0, right=1, data=data_prussia_fin)
summary(servants_ctrl_tobit5)


